/*
  This file is part of code_saturne, a general-purpose CFD tool.

  Copyright (C) 1998-2025 EDF S.A.

  This program is free software; you can redistribute it and/or modify it under
  the terms of the GNU General Public License as published by the Free Software
  Foundation; either version 2 of the License, or (at your option) any later
  version.

  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  details.

  You should have received a copy of the GNU General Public License along with
  this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
  Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/

/*----------------------------------------------------------------------------*/

#include "base/cs_defs.h"

/*----------------------------------------------------------------------------
 * Standard C library headers
 *----------------------------------------------------------------------------*/

#include <assert.h>
#include <errno.h>
#include <ctype.h>
#include <stdio.h>
#include <stdarg.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>

#if defined(HAVE_MPI)
#include <mpi.h>
#endif

/*----------------------------------------------------------------------------
 * Local headers
 *----------------------------------------------------------------------------*/

#include "bft/bft_error.h"
#include "bft/bft_printf.h"

#include "atmo/cs_air_props.h"
#include "base/cs_base.h"
#include "base/cs_boundary_conditions.h"
#include "cdo/cs_domain.h"
#include "base/cs_field.h"
#include "base/cs_field_default.h"
#include "base/cs_field_pointer.h"
#include "base/cs_log.h"
#include "base/cs_math.h"
#include "base/cs_mem.h"
#include "base/cs_scalar_clipping.h"
#include "mesh/cs_mesh.h"
#include "mesh/cs_mesh_location.h"
#include "mesh/cs_mesh_quantities.h"
#include "base/cs_parall.h"
#include "base/cs_parameters.h"
#include "base/cs_parameters_check.h"
#include "base/cs_physical_constants.h"
#include "base/cs_prototypes.h"
#include "base/cs_thermal_model.h"
#include "pprt/cs_physical_model.h"

/*----------------------------------------------------------------------------
 *  Header for the current file
 *----------------------------------------------------------------------------*/

#include "atmo/cs_atmo.h"
#include "atmo/cs_atmo_aerosol_ssh.h"
#include "atmo/cs_atmo_chemistry.h"

/*----------------------------------------------------------------------------*/

BEGIN_C_DECLS

/*! \cond DOXYGEN_SHOULD_SKIP_THIS */

/*=============================================================================
 * Local Macro Definitions
 *============================================================================*/

/*=============================================================================
 * Local Type Definitions
 *============================================================================*/

/* atmo chemistry options structure */
static cs_atmo_chemistry_t _atmo_chem = {
  .model = 0,
  .n_species = 0,
  .n_reactions = 0,
  .chemistry_sep_mode = 2,
  .chemistry_with_photolysis = true,
  .aerosol_model = CS_ATMO_AEROSOL_OFF,
  .frozen_gas_chem = false,
  .init_gas_with_lib = false,
  .init_aero_with_lib = false,
  .n_layer = 0,
  .n_size = 0,
  .spack_file_name = nullptr,
  .species_to_field_id = nullptr,
  .species_profiles_to_field_id = nullptr,
  .molar_mass = nullptr,
  .chempoint = nullptr,
  .conv_factor_jac = nullptr,
  .reacnum = nullptr,
  .dlconc0 = nullptr,
  .aero_file_name = nullptr,
  .chem_conc_file_name = nullptr,
  .aero_conc_file_name = nullptr,
  .nt_step_profiles   = 0,
  .n_z_profiles       = 0,
  .n_species_profiles = 0,
  .conc_profiles   = nullptr,
  .z_conc_profiles = nullptr,
  .t_conc_profiles = nullptr,
  .x_conc_profiles = nullptr,
  .y_conc_profiles = nullptr,
  .dt_chemistry_max = 10.0
};

/*============================================================================
 * Static global variables
 *============================================================================*/

cs_atmo_chemistry_t *cs_glob_atmo_chemistry = &_atmo_chem;

static const char *cs_atmo_aerosol_type_enum_name[]
  = {"CS_ATMO_AEROSOL_OFF",
     "CS_ATMO_AEROSOL_SSH"};

static const char *cs_atmo_aerosol_type_name[]
  = {N_("No atmospheric aerosol"),
     N_("Atmospheric aerosol using external code SSH-aerosol")};

static int _init_atmo_chemistry = 1;

/*============================================================================
 * Prototypes for external Polyphemus functions generated by SPACK.
 *============================================================================*/

void
cs_ext_polyphemus_ssh_dimensions(int  *spack_n_species,
                                 int  *n_reactions,
                                 int  *n_photolysis);

void
cs_ext_polyphemus_fexchem_1(int       nespg,
                            int       nrg,
                            double    dlconc[],
                            double    rk[],
                            double    source[],
                            double    conv_factor[],
                            double    dchema[]);

void
cs_ext_polyphemus_fexchem_2(int       nespg,
                            int       nrg,
                            double    dlconc[],
                            double    rk[],
                            double    source[],
                            double    conv_factor[],
                            double    dchema[]);

void
cs_ext_polyphemus_fexchem_3(int       nespg,
                            int       nrg,
                            double    dlconc[],
                            double    rk[],
                            double    source[],
                            double    conv_factor[],
                            double    dchema[]);

void
cs_ext_polyphemus_fexchem_4(int       nespg,
                            int       nrg,
                            double    dlconc[],
                            double    rk[],
                            double    source[],
                            double    conv_factor[],
                            double    dchema[]);
void
cs_ext_polyphemus_jacdchemdc_1(int        nespg,
                               int        nrg,
                               double     dlconc[],
                               double     conv_factor[],
                               double     conv_factor_jac[],
                               double     dlrki[],
                               double     dldrdc[]);

void
cs_ext_polyphemus_jacdchemdc_2(int        nespg,
                               int        nrg,
                               double     dlconc[],
                               double     conv_factor[],
                               double     conv_factor_jac[],
                               double     dlrki[],
                               double     dldrdc[]);
void
cs_ext_polyphemus_jacdchemdc_3(int        nespg,
                               int        nrg,
                               double     dlconc[],
                               double     conv_factor[],
                               double     conv_factor_jac[],
                               double     dlrki[],
                               double     dldrdc[]);

void
cs_ext_polyphemus_ssh_jacdchemdc(int        nespg,
                                 int        nrg,
                                 double     dlconc[],
                                 double     conv_factor[],
                                 double     conv_factor_jac[],
                                 double     dlrki[],
                                 double     dldrdc[]);

void
cs_ext_polyphemus_lu_decompose_1(int      n_species_g,
                                 double   dlalu[]);

void
cs_ext_polyphemus_lu_solve_1(int       n_species_g,
                             double    dlalu[],
                             double    dlx[]);

void
cs_ext_polyphemus_lu_decompose_2(int      n_species_g,
                                 double   dlalu[]);

void
cs_ext_polyphemus_lu_solve_2(int       n_species_g,
                             double    dlalu[],
                             double    dlx[]);

void
cs_ext_polyphemus_lu_decompose_3(int      n_species_g,
                                 double   dlalu[]);

void
cs_ext_polyphemus_lu_solve_3(int       n_species_g,
                             double    dlalu[],
                             double    dlx[]);

void
cs_ext_polyphemus_ssh_lu_decompose(int       n_species_g,
                                   double    dlalu[]);

void
cs_ext_polyphemus_ssh_lu_solve(int       n_species_g,
                               double    dlalu[],
                               double    dlx[]);

/*============================================================================
 * Private function definitions
 *============================================================================*/

static void
_chem_initialize_species_profiles_to_fid(int *species_profiles_fid)
{
  if (_atmo_chem.species_profiles_to_field_id == nullptr)
    CS_MALLOC(_atmo_chem.species_profiles_to_field_id,
              _atmo_chem.n_species_profiles, int);

  for (int i = 0; i < _atmo_chem.n_species_profiles; i++)
    _atmo_chem.species_profiles_to_field_id[i] = species_profiles_fid[i];
}

/*----------------------------------------------------------------------------*/
/*
 * \brief Beta function: beta(x,y) = gamma(x)*gamma(y)/gamma(x+y)
 */
/*----------------------------------------------------------------------------*/

static inline cs_real_t
_beta(cs_real_t x,
      cs_real_t y)
{
  return tgamma(x)*tgamma(y)/tgamma(x+y);
}

/*----------------------------------------------------------------------------*/
/*
 * \brief Compute hypergeometric function for |x| < 1 using a series
 *        (cf. for the definition of this function, see for example:
 *        http://mathworld.wolfram.com/hypergeometricfunction.html)
 */
/*----------------------------------------------------------------------------*/

static cs_real_t
_hypser(cs_real_t  a,
        cs_real_t  b,
        cs_real_t  c,
        cs_real_t  x)
{
  cs_real_t hypser;

  const int maxiter = 10000;
  const cs_real_t error = 1.e-08;

  if (cs::abs(x) >= 1)
    bft_error(__FILE__, __LINE__, 0,
              "%s\n"
              "The x parameter should verify |x| < 1,  x = %12.5e",
              __func__, x);

  cs_real_t fac = 1;
  cs_real_t temp = fac;
  cs_real_t aa  = a;
  cs_real_t bb  = b;
  cs_real_t cc  = c;

  for (int nn = 1; nn < maxiter; nn++) {
    fac = ((aa*bb)/cc)*fac;
    fac = fac*x/nn;
    hypser = fac + temp;

    if (cs::abs(hypser - temp) <= error)
      return hypser;

    aa += 1;
    bb += 1;
    cc += 1;
    temp = hypser;
  }

  return hypser;
}

/*----------------------------------------------------------------------------*/
/*
 * \brief Hypergeometric function
 *        (see http://mathworld.wolfram.com/hypergeometricfunction.html
 *        for definition)
 */
/*----------------------------------------------------------------------------*/

static cs_real_t
_hypgeo(cs_real_t  a,
        cs_real_t  b,
        cs_real_t  c,
        cs_real_t  x)
{
  cs_real_t hypgeo;

  const cs_real_t pp = 0.1;

  /* Initialization */

  const cs_real_t gammaa   = tgamma(a);
  const cs_real_t gammab   = tgamma(b);
  const cs_real_t gammac   = tgamma(c);
  const cs_real_t gammabma = tgamma(b-a);
  const cs_real_t gammacma = tgamma(c-a);
  const cs_real_t gammaamb = tgamma(a-b);
  const cs_real_t gammacmb = tgamma(c-b);

  /* Compute hypergeometric function by convergent series for |x|<1 */

  if (x >= (pp - 1)) {
    hypgeo = _hypser(a, b, c, x);
                     }
  else if (x >= -1-pp) {
    cs_real_t y1 = _hypser(a, a+1.-c, a+1.-b, 1./x);
    cs_real_t y2 = _hypser(b, b+1.-c, b+1.-a, 1./x);
    hypgeo = (gammac*gammabma*y1*pow(-1./x, a))/(gammab*gammacma)
           + (gammac*gammaamb*y2*pow(-1./x, b))/(gammaa*gammacmb);
  }
  else {
    cs_real_t y1 = _hypser(a, a+1.-c, a+1.-b, 1./(-1.-pp));
    cs_real_t y2 = _hypser(b, b+1.-c, b+1.-a, 1./(-1.-pp));
    cs_real_t hyp1 = (gammac*gammabma*y1*pow(-1./(-1.-pp), a))/(gammab*gammacma)
                   + (gammac*gammaamb*y2*pow(-1./(-1.-pp), b))/(gammaa*gammacmb);
    cs_real_t hyp2 = _hypser(a, b, c, -1.+pp);
    hypgeo = hyp1 + (x - (-1.-pp))*(hyp2 - hyp1)/(2.*pp);
  }

  return hypgeo;
}

/*----------------------------------------------------------------------------
 * Convert string to lower case
 *----------------------------------------------------------------------------*/

static void
_strtolower(char        *dest,
            const char  *src)
{
  char *_dest = dest;
  while (*src) {
    *_dest = tolower(*src);
    src++;
    _dest++;
  }
}

/*----------------------------------------------------------------------------*/
/*
 * \brief  Solver of AX=B with LU decomposition of A for atmospheric chemical
 *         systems
 *
 * \param[in]      kindlu        flag for LU decomposition
 * \param[in]      dla           Matrix A in AX=B
 * \param[in,out]  dlalu         Matrix A in AX=B after LU decomposition
 * \param[in]      dlx           Vector of unknowns X in AX=B
 * \param[in]      dlb           Vector B in AX=B
 */
/*----------------------------------------------------------------------------*/

static void
_decompose_lu(const int  n_species_g,
              const int  kindlu,
              cs_real_t  dla[],
              cs_real_t  dlalu[],
              cs_real_t  dlx[],
              cs_real_t  dlb[])
{
  const cs_atmo_chemistry_t *atmo_chem = cs_glob_atmo_chemistry;

  /* Setup
     ----- */

  for (int ji = 0; ji < n_species_g; ji++)
    dlx[ji] = dlb[ji];

  /* Compute DLx
     ----------- */

  if (kindlu == 0) {
    for (int jj = 0; jj < n_species_g; jj++) {
      for (int ji = 0; ji < n_species_g; ji++) {
        dlalu[ji + jj*n_species_g] = dla[ji + jj*n_species_g];
      }
    }

   if (atmo_chem->model == 1) {
     cs_ext_polyphemus_lu_decompose_1(n_species_g, dlalu);
     cs_ext_polyphemus_lu_solve_1(n_species_g, dlalu, dlx);
   }
   else if (atmo_chem->model == 2) {
     cs_ext_polyphemus_lu_decompose_2(n_species_g, dlalu);
     cs_ext_polyphemus_lu_solve_2(n_species_g, dlalu, dlx);
   }
   else if (atmo_chem->model == 3) {
     cs_ext_polyphemus_lu_decompose_3(n_species_g, dlalu);
     cs_ext_polyphemus_lu_solve_3(n_species_g, dlalu, dlx);
   }
   else if (atmo_chem->model == 4) {
     cs_ext_polyphemus_ssh_lu_decompose(n_species_g, dlalu);
     cs_ext_polyphemus_ssh_lu_solve(n_species_g, dlalu, dlx);
   }

  }
  else {
    if (atmo_chem->model == 1) {
      cs_ext_polyphemus_lu_solve_1(n_species_g, dlalu, dlx);
    }
    else if (atmo_chem->model == 2) {
      cs_ext_polyphemus_lu_solve_2(n_species_g, dlalu, dlx);
    }
    else if (atmo_chem->model == 3) {
      cs_ext_polyphemus_lu_solve_3(n_species_g, dlalu, dlx);
    }
    else if (atmo_chem->model == 4) {
      cs_ext_polyphemus_ssh_lu_solve(n_species_g, dlalu, dlx);
   }

  }
}

/*----------------------------------------------------------------------------*/
/*
 * \brief  Rosenbrock solver for atmospheric chemistry
 *
 * \param[in,out] dlconc        concentrations vector
 * \param[in]     zcsourc       source term for first iteration
 * \param[in]     zcsourcf      source term for second iteration
 * \param[in]     conv_factor   conversion factor
 * \param[in]     dlstep        time step
 * \param[in]     dlrki         kinetic rates for first iteration
 * \param[in]     dlrkf         kinetic rates for second iteration
 */
/*----------------------------------------------------------------------------*/

static void
_rosenbrock_solver(cs_real_t dlconc[],
                   cs_real_t zcsourc[],
                   cs_real_t zcsourcf[],
                   cs_real_t conv_factor[],
                   cs_real_t dlstep,
                   cs_real_t dlrki[],
                   cs_real_t dlrkf[])
{
  const cs_atmo_chemistry_t *atmo_chem = cs_glob_atmo_chemistry;
  const int n_reactions = atmo_chem->n_reactions;
  const int n_species_g = atmo_chem->n_species;

  const cs_real_t igamma = 1.0 + 1.0/sqrt(2.0);

  /* Computes the chemistry
     --------------------- */

  cs_real_t dlr[n_species_g];

  if (atmo_chem->model == 1)
    cs_ext_polyphemus_fexchem_1(n_species_g,
                                n_reactions,
                                dlconc, dlrki,
                                zcsourc, conv_factor, dlr);
  else if (atmo_chem->model == 2)
    cs_ext_polyphemus_fexchem_2(n_species_g,
                                n_reactions,
                                dlconc, dlrki,
                                zcsourc, conv_factor, dlr);
  else if (atmo_chem->model == 3)
    cs_ext_polyphemus_fexchem_3(n_species_g,
                                n_reactions,
                                dlconc, dlrki,
                                zcsourc, conv_factor, dlr);
  else if (atmo_chem->model == 4)
    cs_ext_polyphemus_fexchem_4(n_species_g,
                                n_reactions,
                                dlconc, dlrki,
                                zcsourc, conv_factor, dlr);

  /* Compute the jacobian
     -------------------- */
  cs_real_t dldrdc[n_species_g*n_species_g];

  if (atmo_chem->model == 1)
    cs_ext_polyphemus_jacdchemdc_1(n_species_g,
                                   n_reactions,
                                   dlconc,
                                   conv_factor,
                                   _atmo_chem.conv_factor_jac,
                                   dlrki, dldrdc);
  else if (atmo_chem->model == 2)
    cs_ext_polyphemus_jacdchemdc_2(n_species_g,
                                   n_reactions,
                                   dlconc,
                                   conv_factor,
                                   _atmo_chem.conv_factor_jac,
                                   dlrki, dldrdc);
  else if (atmo_chem->model == 3)
    cs_ext_polyphemus_jacdchemdc_3(n_species_g,
                                   n_reactions,
                                   dlconc,
                                   conv_factor,
                                   _atmo_chem.conv_factor_jac,
                                   dlrki, dldrdc);
  else if (atmo_chem->model == 4)
    cs_ext_polyphemus_ssh_jacdchemdc(n_species_g,
                                     n_reactions,
                                     dlconc,
                                     conv_factor,
                                     _atmo_chem.conv_factor_jac,
                                     dlrki, dldrdc);

   /* Computes K1 system: DLmat * K1 = DLb1
      -------------------------------------- */
   cs_real_t dlb1[n_species_g], dlk1[n_species_g];
   cs_real_t dlmat[n_species_g*n_species_g], dlmatlu[n_species_g*n_species_g];

   for (int jj = 0; jj < n_species_g; jj++) {
     dlb1[jj] = dlr[jj];
     for (int ji = 0; ji < n_species_g; ji++)
       dlmat[ji + jj*n_species_g] = -igamma*dlstep*dldrdc[ji + jj*n_species_g];
     dlmat[jj + jj*n_species_g] = 1.0 + dlmat[jj + jj*n_species_g];
   }

   _decompose_lu(n_species_g, 0, dlmat, dlmatlu, dlk1, dlb1);

   /* Compute K2 system: DLmat * K2 = DLb2
      ------------------------------------- */

   cs_real_t dlconcbis[n_species_g];
   for (int ji = 0; ji < n_species_g; ji++) {
     dlconcbis[ji] = dlconc[ji] + dlstep * dlk1[ji];
     if (dlconcbis[ji] < 0.0) {
       dlconcbis[ji] = 0.0;
       dlk1[ji] = (dlconcbis[ji] - dlconc[ji]) / dlstep;
     }
   }

   if (atmo_chem->model == 1)
     cs_ext_polyphemus_fexchem_1(n_species_g,
                                 n_reactions,
                                 dlconcbis,
                                 dlrkf, zcsourcf,
                                 conv_factor, dlr);
   else if (atmo_chem->model == 2)
     cs_ext_polyphemus_fexchem_2(n_species_g,
                                 n_reactions,
                                 dlconcbis,
                                 dlrkf, zcsourcf,
                                 conv_factor, dlr);
   else if (atmo_chem->model == 3)
     cs_ext_polyphemus_fexchem_3(n_species_g,
                                 n_reactions,
                                 dlconcbis,
                                 dlrkf, zcsourcf,
                                 conv_factor, dlr);
   else if (atmo_chem->model == 4)
     cs_ext_polyphemus_fexchem_4(n_species_g,
                                 n_reactions,
                                 dlconcbis,
                                 dlrkf, zcsourcf,
                                 conv_factor, dlr);

   cs_real_t dlb2[n_species_g], dlk2[n_species_g];
   for (int ji = 0; ji < n_species_g; ji++)
     dlb2[ji] =  dlr[ji] - 2.0*dlk1[ji];

   _decompose_lu(n_species_g, 1, dlmat, dlmatlu, dlk2, dlb2);

   /* Outputs - Compute DLconc - Advance the time
      ------------------------------------------- */

   for (int ji = 0; ji < n_species_g; ji++) {
     dlconc[ji] +=   1.5 * dlstep * dlk1[ji]
                   + 0.5 * dlstep * dlk2[ji];
     if (dlconc[ji] < 0.0)
       dlconc[ji] = 0.0;
   }
}

/*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */

/*============================================================================
 * Public function definitions
 *============================================================================*/

/*----------------------------------------------------------------------------*/
/*!
 * \brief This function initializes the external aerosol code
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_aerosol_initialize(void)
{
  assert(cs_glob_atmo_chemistry->aerosol_model != CS_ATMO_AEROSOL_OFF);

  if (cs_glob_atmo_chemistry->aerosol_model == CS_ATMO_AEROSOL_SSH)
    cs_atmo_aerosol_ssh_initialize();
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief This function finalizes the external aerosol code.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_aerosol_finalize(void)
{
  assert(cs_glob_atmo_chemistry->aerosol_model != CS_ATMO_AEROSOL_OFF);

  if (cs_glob_atmo_chemistry->aerosol_model == CS_ATMO_AEROSOL_SSH)
    cs_atmo_aerosol_ssh_finalize();

  CS_FREE(cs_glob_atmo_chemistry->aero_conc_file_name);
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief This function fills the given array with gas concentrations from
 *        the external aerosol code.
 *
 * \param[out]  array  gas concentrations
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_aerosol_get_gas(cs_real_t  *array)
{
  assert(cs_glob_atmo_chemistry->aerosol_model != CS_ATMO_AEROSOL_OFF);

  if (cs_glob_atmo_chemistry->aerosol_model == CS_ATMO_AEROSOL_SSH)
    cs_atmo_aerosol_ssh_get_gas(array);
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief This function computes a time step of gaseous chemistry and aerosols
 *        dynamic using the external aerosol code.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_aerosol_time_advance(void)
{
  assert(cs_glob_atmo_chemistry->aerosol_model != CS_ATMO_AEROSOL_OFF);

  if (cs_glob_atmo_chemistry->aerosol_model == CS_ATMO_AEROSOL_SSH)
    cs_atmo_aerosol_ssh_time_advance();
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief Compute aerosol cloud droplets nucleation when using the atmospheric
 * humid model using a microphysical model.
 *
 * It is taken into account as an additional step split from advection-diffusion
 * equation, hence the droplet number is first clipped if necessary.
 *
 * \param[out]  nc      droplet number (scalar) in 1/cm**3
 * \param[in]   rom     density of air in kg/m**3
 * \param[in]   qldia   mass fraction of liquid water in kg/kg
 * \param[in]   pphy    true pressure in pascals
 * \param[in]   refrad  radiative cooling
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_aerosol_nuclea(cs_real_t         *nc,
                       const cs_real_t   *rom,
                       const cs_real_t   *qldia,
                       const cs_real_t   *pphy,
                       const cs_real_t   *refrad)
{
  /* Initialization
   * -------------- */

  const cs_mesh_t *mesh = cs_glob_mesh;
  const cs_lnum_t n_cells = mesh->n_cells;

  const cs_real_3_t *vel = (const cs_real_3_t *)CS_F_(vel)->val;
  const cs_real_t *tempc = cs_field_by_name("real_temperature")->val;

  cs_fluid_properties_t *phys_pro = cs_get_glob_fluid_properties();

  cs_real_t cp = phys_pro->cp0;
  cs_real_t rvsra = phys_pro->rvsra;
  cs_real_t clatev = phys_pro->clatev;
  cs_real_t rair = phys_pro->r_pg_cnst;

  const cs_real_t rhowater = 1000; // kg/m^3
  const cs_real_t tkelvi = cs_physical_constants_celsius_to_kelvin;

  cs_real_t nuc = 0.;
  cs_real_t fbeta = 0.;
  cs_real_t constc = 0.;
  cs_real_t constk = 0.;
  cs_real_t constmu = 0;
  cs_real_t constbeta = 0;
  cs_real_t rvap = rair*rvsra;
  const int modnuc = cs_glob_atmo_option->nucleation_model;

  if (modnuc == 1) {
    /* Constants for the model of Pruppacher and Klett (1997)
     * Case of Buffalo, New-York, Kocmond [1965] */
    constc = 3500.;
    constk = 0.9;

    /* fbeta : beta function (constk/2, 3/2) */
    fbeta = _beta(constk/2., 1.5);
  }
  else if (modnuc == 2) {
    /* Constants for model of Cohard and Pinty (1998)
     * (general case) */
    constc    = 3270.;
    constk    = 1.56;
    constmu   = 0.7;
    constbeta = 136.;

    /* Polluted Air Case. See Hudson and Li (1995)
     * constc    = 1865. ! 8700.d0 if Cohard and Pinty used on PARISFOG
     * constk    = 0.86
     * constmu   = 1.50
     * constbeta = 6.80 */
    fbeta = _beta(constk/2., 1.5);
  }

  /* Computation new Nc field (in cm^-3)
   * ----------------------------------- */

  for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {

    if (qldia[c_id] <= cs_math_epzero) {

      nc[c_id] = 0.;

    }
    else { /* qldia[c_id] > cs_math_epzero */

      /* New droplets are created by nucleation if vel_z > 0 or refrad < 0,
       * then if number of new droplets is superior to nc,
       * the difference is added */

      if ((vel[c_id][2] > cs_math_epzero) || (refrad[c_id] < cs_math_epzero)) {

        cs_real_t tempk = tempc[c_id] + tkelvi;
        cs_real_t esat = cs_air_pwv_sat(tempc[c_id]);
        cs_real_t kka = ((5.69 + 0.017*tempc[c_id])/0.239)*1.e-3;
        cs_real_t ddv = 0.211*pow(tempk/tkelvi, 1.94)*(101325./pphy[c_id])*1.e-4;
        cs_real_t aa1 =   0.622*clatev*9.81/(rair*cp*cs_math_pow2(tempk))
                        - 9.81/(rair*tempk);
        cs_real_t aa2 =   rair*tempk/(0.622*esat)
                        + (0.622*cs_math_pow2(clatev))/(tempk*pphy[c_id]*cp);
        cs_real_t aa3 = 1.0/((1000.*rvap*tempk)/(esat*ddv)
                        + clatev*1000.*(clatev/(tempk*rvap) -1.0)/(kka*tempk));
        cs_real_t aa4 = -0.622*clatev/(rair*cs_math_pow2(tempk));

        if ((aa1*vel[c_id][2]+aa4*refrad[c_id]) > cs_math_epzero) {

          // 1.1  Model of Pruppacher and Klett (1997)
          if (modnuc == 1) {
            const cs_real_t cons2 = pow(constc, 2./(constk+2.));
            const cs_real_t cons3 = pow(0.01*(  aa1*vel[c_id][2]
                                              + aa4*refrad[c_id]), 3./2.);
            const cs_real_t cons4
              = (2.*cs_math_pi*rhowater*aa2*pow(aa3, 3./2.)*constk*fbeta);
            nuc = cons2*pow(cons3/cons4, constk/(constk + 2.));
          }

          // 1.2  Modele of Cohard and Pinty (1998)
          else if (modnuc == 2) {
            /* compute max supersaturation: sursat */
            cs_real_t yy = 1.0;
            cs_real_t sursat = 0;
            cs_real_t tmpsur = 0;
            for (int ii = 0; ii < 20; ii++) {
              tmpsur = sursat;
              yy = _hypgeo(constmu, constk/2, constk/2 + 3./2,
                        -constbeta*cs_math_pow2(sursat));
              const cs_real_t cons1
                = pow(0.01*(aa1*vel[c_id][2] + aa4*refrad[c_id]), 3./2);
              const cs_real_t cons2
                = (2*constk*constc*cs_math_pi*rhowater*aa2*fbeta
                   *pow(aa3, 3.0/2));

              sursat = pow((cons1/cons2)/yy, 1./(constk + 2.));
            }

            if (cs::abs(tmpsur - sursat) > 1.e-2)
              bft_error(__FILE__, __LINE__, 0,
                        "Maximum saturation has not converged\n"
                        "Residue = %12.5e", cs::abs(tmpsur - sursat));

            nuc = constc*pow(sursat, constk)*_hypgeo(constmu,
                                                     constk/2,
                                                     constk/2 + 1.,
                                                     -constbeta*
                                                     cs_math_pow2(sursat));
            if (nuc < 0.)
              bft_error(__FILE__, __LINE__, 0,
                        _("Cohard and Pindy model (1998).\n"
                          "The nucleation is negative."));
          }

          /* 1.3  Model of Abdul-Razzak and al. (1998)
           * This model requires a fine knowledge of aerosols chemistry. */
          else if (modnuc == 3) {

            /* Warning: this set of constants fits for PARISFOG case.
             * They were determined using fine data on aerosols measured
             * during PARISFOG POI 13 */
            cs_real_t sigaero3 = 1.3;
            cs_real_t sigaero2 = 1.59;
            cs_real_t sigaero1 = 1.570;
            cs_real_t ntotal1 = 8700.;
            cs_real_t ntotal2 = 8300.;
            cs_real_t ntotal3 = 1000.;
            cs_real_t rayonm3 = 0.4e-6;
            cs_real_t rayonm2 = 0.055e-6;
            cs_real_t rayonm1 = 0.0165e-6;

            /* surface tension water-steam [N/m], tempk [K] */
            const cs_real_t tauvw = 0.0750;

            /* fractions in %
             * fisrt component: sulfate
             * second component: nitrate
             * third component: black carbon
             * fourth component: organic carbon */
            const cs_real_t fraero[4] = {0.25, 0.2, 0.164, 0.386};

            /* number of ions produced by dissociation of a salt molecule
               in water (if NaCl, nbion=2) */
            const cs_real_t nbion[4] = {3., 2., 0., 1.};

            /* osmotic coefficient */
            const cs_real_t coefosm[4] = {1., 1., 1., 1.};

            /* mass fraction of of soluble
               substance in aerosol mix: msolu/mmix */
            const cs_real_t frasolu[4] = {1., 1., 0., 0.1};

            /* molecular mass of water */
            const cs_real_t mmh2o = 18.e-3;

            /* molecular mass of aerosol components */
            const cs_real_t mmaero[4] = {132e-3, 80e-3, 250e-3, 250e-3};

            /* density of aerosol !FIXME to check */
            const cs_real_t rhoaero[4] = {1.77e3, 1.77e3, 1.77e3, 1.77e3};

            /* coefficient for Kelvin effect: coefa [/m] */
            const cs_real_t coefa = 2*tauvw/(rhowater*rvap*tempk);

            /* FIXME
             * cf Pruppacher and Klett 1997 (reprinted correction 2000) (6-28)
             * coefa = 3.3d-7 / tempk */

            /* coefficient for Raoult effect: coefb [-] */
            cs_real_t numcb = 0;
            cs_real_t dencb = 0;
            for (int ii = 0; ii < 4; ii++) {
              numcb += fraero[ii]*nbion[ii]*coefosm[ii]*frasolu[ii]/mmaero[ii];
              dencb += fraero[ii]/rhoaero[ii];
            }
            const cs_real_t coefb = mmh2o*numcb/(dencb*rhowater);

            /* supersaturation [-] */
            cs_real_t sursat1 = (2./sqrt(coefb))*(pow(coefa/3./rayonm1, 1.5));
            cs_real_t sursat2 = (2./sqrt(coefb))*(pow(coefa/3./rayonm2, 1.5));
            cs_real_t sursat3 = (2./sqrt(coefb))*(pow(coefa/3./rayonm3, 1.5));

            /* standard deviation function */
            cs_real_t foncf1 = 0.5*exp(2.5*cs_math_pow2(log(sigaero1)));
            cs_real_t foncf2 = 0.5*exp(2.5*cs_math_pow2(log(sigaero2)));
            cs_real_t foncf3 = 0.5*exp(2.5*cs_math_pow2(log(sigaero3)));
            cs_real_t foncg1 = 1.0 + 0.25*(log(sigaero1));
            cs_real_t foncg2 = 1.0 + 0.25*(log(sigaero2));
            cs_real_t foncg3 = 1.0 + 0.25*(log(sigaero3));

            /* coefficient corresponding to vertical velocity |aa1| */
            aa1 =   0.622*clatev*9.81/(rair*cp*cs_math_pow2(tempk))
                  - 9.81/(rair*tempk);

            /* coefficient corresponding to liquid water condensation |aa2| */
            aa2 =   rair*tempk/(0.622*esat)+(0.622*cs_math_pow2(clatev))
                  / (tempk*pphy[c_id]*cp);

            /* coefficient corresponding to the droplet growth |aa3| */
            ddv = 0.211*pow(tempk/tkelvi, 1.94)*(101325./pphy[c_id])*1.e-4;
            rvap = rvsra*rair;
            aa3 =   1/((1000*rvap*tempk)/(esat*ddv)
                  + clatev*1000*(clatev/(tempk*rvap)-1.0)/(kka*tempk));

            /* coefficient corresponding to infrared radiation |aa4| */
            aa4 = -0.622*clatev/(rair*cs_math_pow2(tempk));

            /* alphaV/G */
            cs_real_t coefavg = (aa1*vel[c_id][2] + aa4*refrad[c_id])/aa3;

            /* coefficient zeta */
            cs_real_t coefzeta = 2*coefa/3*sqrt(coefavg);

            /* coefficient eta - Ntot [m^(-3)] */
            const cs_real_t coefeta1
              = pow(coefavg, 1.5)/(2*cs_math_pi*rhowater*aa2*ntotal1*1.e6);
            const cs_real_t coefeta2
              = pow(coefavg, 1.5)/(2*cs_math_pi*rhowater*aa2*ntotal2*1.e6);
            const cs_real_t coefeta3
              = pow(coefavg, 1.5)/(2*cs_math_pi*rhowater*aa2*ntotal3*1.e6);

            /* max supersaturation (smax) */
            cs_real_t x_tp1 = foncf1*pow(coefzeta/coefeta1, 1.5);
            cs_real_t x_tp2 = pow((sursat1*sursat1)/(coefeta1+3*coefzeta), 0.75);
            cs_real_t smax1 = (x_tp1 + foncg1*x_tp2)/sursat1/sursat1;

            x_tp1 = foncf2*pow(coefzeta/coefeta2, 1.5);
            x_tp2 = pow((sursat2*sursat2)/(coefeta2+3*coefzeta), 0.75);
            cs_real_t smax2 = (x_tp1 + foncg2*x_tp2)/sursat2/sursat2;

            x_tp1 = foncf3*(pow(coefzeta/coefeta3, 1.5));
            x_tp2 = pow((sursat3*sursat3)/(coefeta3+3*coefzeta), 0.75);
            cs_real_t smax3 = (x_tp1 + foncg3*x_tp2)/sursat3/sursat3;

            cs_real_t smax = 1.0/sqrt(smax3 + smax2 + smax1);

            if ((smax > 1.0) || (sursat1 > 1.0) ||
                (sursat2 > 1.0) || (sursat3 > 1.0)) {
              bft_error(__FILE__, __LINE__, 0,
                        _("%s: Abdul-Razzak and al. model (1998).\n"
                          " Negative supersaturation."), __func__);
            }

            cs_real_t f3_ov_sqrt2 = 3.*sqrt(2.);
            cs_real_t ugauss1 =   (2.*log(sursat1/smax))
                                / (f3_ov_sqrt2*log(sigaero1));
            cs_real_t ugauss2 =   (2.*log(sursat2/smax))
                                / (f3_ov_sqrt2*log(sigaero2));
            cs_real_t ugauss3 =   (2.*log(sursat3/smax))
                                / (f3_ov_sqrt2*log(sigaero3));

            const cs_real_t nuc1 = 0.5 * ntotal1 * (1.0-erf(ugauss1));
            const cs_real_t nuc2 = 0.5 * ntotal2 * (1.0-erf(ugauss2));
            const cs_real_t nuc3 = 0.5 * ntotal3 * (1.0-erf(ugauss3));

            nuc = nuc1 + nuc2 + nuc3;

          } /* end nucleation model */

          /* 2. Compute difference */
          nuc = cs::max(nuc - nc[c_id], 0.0);
        }
        else {
          nuc = 0.0;
        }

        nc[c_id] += nuc;
      }

      /* 3. if qc > 0, w <= 0 and nc = 0,
       * we impose nc > 0 so that the mean volume radius = 10 microns */
      else if (nc[c_id] < cs_math_epzero) {

        nc[c_id] =   1.e-6*(3.*rom[c_id]*qldia[c_id])
                   / (4.*cs_math_pi*rhowater*cs_math_pow3(10e-6));

      } /* end Vel_z > cs_math_epzero*/

    } /* end qldia[c_id] > cs_math_epzero */
  }
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief Reads initial aerosol concentration and number
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_read_aerosol(void)
{
  const cs_atmo_chemistry_t *at_chem = cs_glob_atmo_chemistry;
  const cs_fluid_properties_t *phys_pro = cs_get_glob_fluid_properties();

  const int n_aer = at_chem->n_size;
  const int nlayer_aer = at_chem->n_layer;
  const int size = n_aer * nlayer_aer;

  cs_atmo_chem_initialize_dlconc0();

  bft_printf("reading of aerosols numbers and concentrations\n");

  if (at_chem->init_aero_with_lib) {

    // The external library provides the concentrations / numbers
    if (at_chem->aerosol_model == CS_ATMO_AEROSOL_SSH)
      cs_atmo_aerosol_ssh_get_aero(at_chem->dlconc0);

    // Conversion from molecules / m^3 to ppm
    const cs_real_t ro0 = 1e3*phys_pro->ro0;
    for (int jb = 0; jb < size; jb++)
      at_chem->dlconc0[jb] = at_chem->dlconc0[jb] / ro0;

    // Conversion from molecules / m^3 to molecules / kg
    const int end_size = size + n_aer;
    for (int jb = size; jb < end_size; jb++)
      at_chem->dlconc0[jb] = at_chem->dlconc0[jb] / phys_pro->ro0;

  }
  else {

    // Read from file
    FILE *file = fopen(at_chem->aero_conc_file_name, "r");

    // Reading aerosol numbers
    for (int jb = 0; jb < n_aer; jb++)
      fscanf(file,"%le", &at_chem->dlconc0[jb + size]);

    // Reading aerosol concentrations
    for (int jb = 0; jb < n_aer; jb++)
      for (int jsp = 0; jsp < nlayer_aer; jsp++)
        fscanf(file,"%le", &at_chem->dlconc0[jb + jsp*n_aer]);

    fclose(file);

  }

  /* Logging
     ------- */

  bft_printf("\n===================================================\n");
  bft_printf("printing aerosol numbers\n");
  for (int jb = 0; jb < n_aer; jb++) {
    const int f_id
      = at_chem->species_to_field_id[at_chem->n_species + size + jb];
    const cs_field_t *f = cs_field_by_id(f_id);
    bft_printf("%s : %10.2le\n ",
               cs_field_get_label(f),
               at_chem->dlconc0[jb + size]);
  }

  bft_printf("\n===================================================\n");
  bft_printf("printing aerosol concentrations\n");
  for (int jb = 0; jb < n_aer; jb++) {
    bft_printf("Size bin id %d\n", jb);
    for (int jsp = 0; jsp < nlayer_aer; jsp++) {
      const int f_id
        = at_chem->species_to_field_id[at_chem->n_species + jsp*n_aer];
      const cs_field_t *f = cs_field_by_id(f_id);
      bft_printf("%s : %10.2le\n ",
                 cs_field_get_label(f),
                 at_chem->dlconc0[jb + jsp*n_aer]);
    }
  }
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief Reads the chemistry profile data for the atmospheric chemistry
 *
 * \param[in] mode    if false reading of dimensions only else reading of data
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_read_chemistry_profile(int mode)
{
  cs_atmo_option_t *at_opt = cs_glob_atmo_option;
  const int n_fields = cs_field_n_fields();

  FILE *file = nullptr;
  char chaine[100] = "";
  static bool switch_to_labels;
  int species_profiles_to_fid[_atmo_chem.n_species_profiles];
  const char *name = _atmo_chem.chem_conc_file_name;

  if (name != nullptr)
    file = fopen(name, "r");

  if (mode == 1) {

    bft_printf("Reading concentration profiles data.\n");

    for (int itp = 0; itp < _atmo_chem.nt_step_profiles; itp++) {

      int year, quant, hour, minute, second;

      // read time
      while (strcmp(chaine, "second") != 0)
        fscanf(file, "%s", chaine);
      fscanf(file, "%d %d %d %d %d",
             &year, &quant, &hour, &minute, &second);

      if (second < 0 || quant > 366)
        cs_parameters_error
          (CS_ABORT_IMMEDIATE,
           _("Error opening the chemistry profile file"),
           _("verify input format (integer, real)\n"
             "The computation will not be run"));

      if (at_opt->syear < 0.0) {
        at_opt->syear  = year;
        at_opt->squant = quant;
        at_opt->shour  = hour;
        at_opt->smin   = minute;
        at_opt->ssec   = (cs_real_t)second;
      }

      /* Compute the julian day for the starting day of the simulation
       * (julian day at 12h) */
      const cs_real_t sjday
        = at_opt->squant + ( ( 1461 * (at_opt->syear + 4800 + (1 - 14) / 12))/4
                             + (367 * (1 - 2 - 12 * ((1 - 14) / 12))) / 12
                             - (3 * ((at_opt->syear + 4900 + (1-14)/12) / 100))/4
                             + 1 - 32075 ) - 1;
      const cs_real_t jday
        = quant + (  (1461*(year+4800+(1-14)/12))/4
                     + (367*(1-2-12*((1-14) / 12))) / 12
                     - (3 * ((year+4900+(1-14)/12)/100)) / 4
                     +  1-32075) - 1;

      _atmo_chem.t_conc_profiles[itp] = (jday - sjday)*86400.0
        + (hour - at_opt->shour)*3600.0
        + (minute - at_opt->smin)*60.0  + (second - at_opt->ssec);

      if (itp > 0)
        if (_atmo_chem.t_conc_profiles[itp] < _atmo_chem.t_conc_profiles[itp-1])
                  cs_parameters_error
                    (CS_ABORT_IMMEDIATE,
                     _("Error in the chemistry profile file:\n"),
                     _("check that the chronogical order of"
                       "the profiles are respected"));

      while (strcmp(chaine, "domaine") != 0)
        fscanf(file, "%s", chaine);

      /* reading the position of the profile
         -----------------------------------*/
      fscanf(file, "%lf %lf",
             &_atmo_chem.x_conc_profiles[itp],
             &_atmo_chem.y_conc_profiles[itp]);

      char labels[100] = "";

      if (_atmo_chem.init_gas_with_lib) {

        /* read the concentrations
           ----------------------- */
        const cs_fluid_properties_t *phys_pro = cs_get_glob_fluid_properties();
        const cs_real_t ro0 = (1e3*phys_pro->ro0);

        _atmo_chem.z_conc_profiles[0] = 0.0;
        cs_atmo_aerosol_get_gas(_atmo_chem.conc_profiles);

        for (int ii = 0; ii < _atmo_chem.n_species_profiles; ii++)
          _atmo_chem.conc_profiles[ii] = _atmo_chem.conc_profiles[ii]/ro0;

      }
      else {
        if (_atmo_chem.n_species_profiles > 0) {
          if (switch_to_labels) {

            while (strcmp(chaine, "initialisées") != 0)
              fscanf(file, "%s", chaine);
            int nbchim = 0;
            fscanf(file, "%d", &nbchim);
            fscanf(file, "%s", chaine);
            while (strcmp(chaine, "initialisées") != 0)
              fscanf(file, "%s", chaine);
            for (int ii = 0; ii < _atmo_chem.n_species_profiles; ii++) {
              fscanf(file, "%s", labels);
              if (strcmp(labels, "") == 0) {
                cs_parameters_error
                  (CS_ABORT_IMMEDIATE,
                   _("ATMOSPHERIC CHEMISTRY FROM SPACK:\n"),
                   _("Could not identify the given species label\n"
                     "Given species label :%s"), labels);
              }
              else {
                for (int f_id = 0; f_id < n_fields; f_id++) {
                  const cs_field_t *f = cs_field_by_id(f_id);
                  const char *f_label = cs_field_get_label(f);
                  if (strcmp(labels, f_label) != 0)
                    continue;
                  species_profiles_to_fid[ii] = f_id;
                }
              }
            }
            _chem_initialize_species_profiles_to_fid
              (species_profiles_to_fid);
          }

          int nbchmz = 0;
          cs_real_t zconctemp[_atmo_chem.n_species_profiles+1];

          while (strcmp(chaine, "concentrations") != 0)
            fscanf(file, "%s", chaine);
          fscanf(file, "%d", &nbchmz);

          nbchmz = _atmo_chem.n_z_profiles;
          const int size = nbchmz*_atmo_chem.nt_step_profiles;
          for (int ii = 0; ii < _atmo_chem.n_z_profiles; ii++) {
            for (int jj = 0; jj < _atmo_chem.n_species_profiles+1; jj++)
              fscanf(file, "%lf", &zconctemp[jj]);
            _atmo_chem.z_conc_profiles[ii] = zconctemp[0];

            for (int kk = 1; kk < _atmo_chem.n_species_profiles+1; kk++)
              _atmo_chem.conc_profiles[ii+(itp)*nbchmz+(kk-1)*size]
                = zconctemp[kk];
          }
        } // fin test nespgi

      } //fin test init_gas_with_lib

      /* logging and initialaze species profiles
         --------------------------------------- */

      if (itp == 0) {
        bft_printf("===================================================\n");
        bft_printf("Concentration profiles\n");
      }

      bft_printf("year, quant-day, hour, minute, second\n");
      bft_printf("%d  %d  %d  %d  %d\n", year, quant, hour, minute, second);
      bft_printf("t_conc_profiles: %10.4lf\n", _atmo_chem.t_conc_profiles[itp]);

      bft_printf("\n");
      const int nbchmz = _atmo_chem.n_z_profiles;
      const int size = nbchmz*_atmo_chem.nt_step_profiles;
      bft_printf("zproc,");
      for (int kk = 0; kk < _atmo_chem.n_species_profiles; kk++) {
        cs_field_t *f = cs_field_by_id(species_profiles_to_fid[kk]);
        const char *f_label = cs_field_get_label(f);
        bft_printf(" %s,",f_label);
      }
      bft_printf("\n");
      for (int ii = 0; ii < _atmo_chem.n_z_profiles; ii++) {
        bft_printf("%10.2lf ",  _atmo_chem.z_conc_profiles[ii]);
        for (int kk = 0; kk < _atmo_chem.n_species_profiles; kk++) {
          bft_printf("%10.5lf ",
                     _atmo_chem.conc_profiles[ii+(itp)*nbchmz+(kk)*size]);

        }
        bft_printf("\n");
      }

    } // end loop on time

  }
  else {

    _atmo_chem.nt_step_profiles = 1;
    bft_printf("Reading dimensions for concentration profiles.\n");

    if (_atmo_chem.init_gas_with_lib) {

      _atmo_chem.n_z_profiles = 1;
      _atmo_chem.n_species_profiles = _atmo_chem.n_species;
      return;

    }
    else {
      while (strcmp(chaine, "initialisées") != 0)
        fscanf(file, "%s", chaine);
      fscanf(file, "%d", &_atmo_chem.n_species_profiles);

      if (_atmo_chem.n_species_profiles < 0) {
        switch_to_labels = true;
        _atmo_chem.n_species_profiles = -_atmo_chem.n_species_profiles;
      }
      else
        switch_to_labels = false;

      while (strcmp(chaine, "concentrations") != 0)
        fscanf(file, "%s", chaine);
      fscanf(file, "%d", &_atmo_chem.n_z_profiles);
      if (_atmo_chem.n_z_profiles < 2)
        cs_parameters_error
                    (CS_ABORT_IMMEDIATE,
                     _("Error in the chemistry profile file:\n"),
                     _("the number of concentrations measurements\n"
                       " must be larger than 2"));
    }

  }

  if (file != nullptr)
    fclose(file);
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief This function set the file name of the chemistry file.
 *
 * \param[in] file_name  name of the file.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_set_chem_conc_file_name(const char  *file_name)
{
  if (file_name == nullptr) {
    return;
  }

  if (_atmo_chem.chem_conc_file_name == nullptr) {
    CS_MALLOC(_atmo_chem.chem_conc_file_name,
              strlen(file_name) + 1,
              char);
  }
  else {
    CS_REALLOC(_atmo_chem.chem_conc_file_name,
               strlen(file_name) + 1,
               char);
  }

  strcpy(_atmo_chem.chem_conc_file_name, file_name);
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief This function set the file name of the aerosol file.
 *
 * \param[in] file_name  name of the file.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_set_aero_conc_file_name(const char  *file_name)
{
  if (file_name == nullptr) {
    return;
  }
  if (cs_glob_atmo_chemistry->aerosol_model == CS_ATMO_AEROSOL_OFF) {
    return;
  }

  if (_atmo_chem.aero_conc_file_name == nullptr) {
    CS_MALLOC(_atmo_chem.aero_conc_file_name,
              strlen(file_name) + 1,
              char);
  }
  else {
    CS_REALLOC(_atmo_chem.aero_conc_file_name,
               strlen(file_name) + 1,
               char);
  }

  strcpy(_atmo_chem.aero_conc_file_name, file_name);
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief Initialize chemistry array.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_init_chemistry(void)
{
  /* Initialization of the chemical scheme
     quasi steady equilibrium NOx scheme */
  if (_atmo_chem.model == 1) {

    _atmo_chem.n_species = 4;
    _atmo_chem.n_reactions = 5;

    if (_atmo_chem.chempoint == nullptr)
      CS_MALLOC(_atmo_chem.chempoint, _atmo_chem.n_species, int);

    if (_atmo_chem.molar_mass == nullptr)
      CS_MALLOC(_atmo_chem.molar_mass, _atmo_chem.n_species, cs_real_t);

    _atmo_chem.chempoint[0] = 4;
    _atmo_chem.chempoint[1] = 3;
    _atmo_chem.chempoint[2] = 2;
    _atmo_chem.chempoint[3] = 1;
    _atmo_chem.molar_mass[0] = 30.0;
    _atmo_chem.molar_mass[1] = 46.0;
    _atmo_chem.molar_mass[2] = 48.0;
    _atmo_chem.molar_mass[3] = 16.0;

    if (_atmo_chem.species_to_field_id == nullptr)
      CS_MALLOC(_atmo_chem.species_to_field_id, _atmo_chem.n_species, int);

    const char *label[] = {"NO", "NO2", "O3", "O3P"};

    const char *name[]
      = {"species_no", "species_no2", "species_o3", "species_o3p"};

    for (int ii = 0; ii < _atmo_chem.n_species; ii++) {
      int f_id = cs_variable_field_create(name[ii],
                                          label[ii],
                                          CS_MESH_LOCATION_CELLS,
                                          1);
      cs_field_t *f = cs_field_by_id(f_id);
      cs_add_model_field_indexes(f);
      _atmo_chem.species_to_field_id[ii] = f->id;
    }

  }
  else if (_atmo_chem.model == 2) {

    _atmo_chem.n_species = 20;
    _atmo_chem.n_reactions = 34;

    if (_atmo_chem.chempoint == nullptr)
      CS_MALLOC(_atmo_chem.chempoint, _atmo_chem.n_species, int);
    if (_atmo_chem.molar_mass == nullptr)
      CS_MALLOC(_atmo_chem.molar_mass, _atmo_chem.n_species, cs_real_t);

    _atmo_chem.chempoint[0]  = 20, _atmo_chem.chempoint[1]  = 19;
    _atmo_chem.chempoint[2]  = 16, _atmo_chem.chempoint[3]  = 17;
    _atmo_chem.chempoint[4]  = 2,  _atmo_chem.chempoint[5]  = 15;
    _atmo_chem.chempoint[6]  = 14, _atmo_chem.chempoint[7]  = 3;
    _atmo_chem.chempoint[8]  = 18, _atmo_chem.chempoint[9]  = 7;
    _atmo_chem.chempoint[10] = 8,  _atmo_chem.chempoint[11] = 9;
    _atmo_chem.chempoint[12] = 4,  _atmo_chem.chempoint[13] = 10;
    _atmo_chem.chempoint[14] = 1,  _atmo_chem.chempoint[15] = 12;
    _atmo_chem.chempoint[16] = 11, _atmo_chem.chempoint[17] = 13;
    _atmo_chem.chempoint[18] = 5,  _atmo_chem.chempoint[19] = 6;

    _atmo_chem.molar_mass[0]  = 30.000;  // Molar mass (g/mol) NO
    _atmo_chem.molar_mass[1]  = 46.000;  // Molar mass (g/mol) NO2
    _atmo_chem.molar_mass[2]  = 48.000;  // Molar mass (g/mol) O3
    _atmo_chem.molar_mass[3]  = 16.000;  // Molar mass (g/mol) O3P
    _atmo_chem.molar_mass[4]  = 16.000;  // Molar mass (g/mol) O1D
    _atmo_chem.molar_mass[5]  = 17.010;  // Molar mass (g/mol) OH
    _atmo_chem.molar_mass[6]  = 33.010;  // Molar mass (g/mol) HO2
    _atmo_chem.molar_mass[7]  = 34.010;  // Molar mass (g/mol) H2O2
    _atmo_chem.molar_mass[8]  = 62.010;  // Molar mass (g/mol) NO3
    _atmo_chem.molar_mass[9]  = 108.01;  // Molar mass (g/mol) N2O5
    _atmo_chem.molar_mass[10] = 47.010;  // Molar mass (g/mol) HONO
    _atmo_chem.molar_mass[11] = 63.010;  // Molar mass (g/mol) HNO3
    _atmo_chem.molar_mass[12] = 28.010;  // Molar mass (g/mol) CO
    _atmo_chem.molar_mass[13] = 30.030;  // Molar mass (g/mol) HCHO
    _atmo_chem.molar_mass[14] = 44.050;  // Molar mass (g/mol) ALD2
    _atmo_chem.molar_mass[15] = 75.040;  // Molar mass (g/mol) C2O3
    _atmo_chem.molar_mass[16] = 121.05;  // Molar mass (g/mol) PAN
    _atmo_chem.molar_mass[17] = 47.030;  // Molar mass (g/mol) XO2
    _atmo_chem.molar_mass[18] = 64.060;  // Molar mass (g/mol) SO2
    _atmo_chem.molar_mass[19] = 98.080;  // Molar mass (g/mol) H2SO4

    const char *label[] = {"NO",   "NO2",  "O3",   "O3P",  "O1D",
                           "OH",   "HO2",  "H2O2", "NO3",  "N2O5",
                           "HONO", "HNO3", "CO",   "HCHO", "ALD2",
                           "C2O3", "PAN",  "XO2",  "SO2",  "H2SO4"};

    const char *name[]
      = {"species_no",  "species_no2",  "species_o3",  "species_o3p",
         "species_o1d", "species_oh",   "species_ho2", "species_h2o2",
         "species_no3", "species_n2o5", "species_hono", "species_hno3",
         "species_co",  "species_hcho", "species_ald2","species_c2o3",
         "species_pan", "species_xo2",  "species_so2", "species_h2so4"};

    if (_atmo_chem.species_to_field_id == nullptr)
      CS_MALLOC(_atmo_chem.species_to_field_id, _atmo_chem.n_species, int);

    for (int ii = 0; ii < _atmo_chem.n_species; ii++) {
      int f_id = cs_variable_field_create(name[ii],
                                          label[ii],
                                          CS_MESH_LOCATION_CELLS,
                                          1);
      cs_field_t *f = cs_field_by_id(f_id);
      cs_add_model_field_indexes(f);
      _atmo_chem.species_to_field_id[ii] = f->id;
    }

  }
  else if (_atmo_chem.model == 3) {
    _atmo_chem.n_species = 52;
    _atmo_chem.n_reactions = 155;

    if (_atmo_chem.chempoint == nullptr)
      CS_MALLOC(_atmo_chem.chempoint, _atmo_chem.n_species, int);
    if (_atmo_chem.molar_mass == nullptr)
      CS_MALLOC(_atmo_chem.molar_mass, _atmo_chem.n_species, cs_real_t);

    _atmo_chem.chempoint[0]  = 48, _atmo_chem.chempoint[1]  = 52;
    _atmo_chem.chempoint[2]  = 47, _atmo_chem.chempoint[3]  = 43;
    _atmo_chem.chempoint[4]  = 1,  _atmo_chem.chempoint[5]  = 42;
    _atmo_chem.chempoint[6]  = 50, _atmo_chem.chempoint[7]  = 17;
    _atmo_chem.chempoint[8]  = 44, _atmo_chem.chempoint[9]  = 9;
    _atmo_chem.chempoint[10] = 15, _atmo_chem.chempoint[11] = 38;
    _atmo_chem.chempoint[12] = 13, _atmo_chem.chempoint[13] = 37;
    _atmo_chem.chempoint[14] = 41, _atmo_chem.chempoint[15] = 45;
    _atmo_chem.chempoint[16] = 51, _atmo_chem.chempoint[17] = 10;
    _atmo_chem.chempoint[18] = 35, _atmo_chem.chempoint[19] = 46;
    _atmo_chem.chempoint[20] = 14, _atmo_chem.chempoint[21] = 49;
    _atmo_chem.chempoint[22] = 39, _atmo_chem.chempoint[23] = 33;
    _atmo_chem.chempoint[24] = 2,  _atmo_chem.chempoint[25] = 3;
    _atmo_chem.chempoint[26] = 40, _atmo_chem.chempoint[27] = 11;
    _atmo_chem.chempoint[28] = 19, _atmo_chem.chempoint[29] = 20;
    _atmo_chem.chempoint[30] = 4,  _atmo_chem.chempoint[31] = 21;
    _atmo_chem.chempoint[32] = 36, _atmo_chem.chempoint[33] = 22;
    _atmo_chem.chempoint[34] = 34, _atmo_chem.chempoint[35] = 16;
    _atmo_chem.chempoint[36] = 23, _atmo_chem.chempoint[37] = 24;
    _atmo_chem.chempoint[38] = 25, _atmo_chem.chempoint[39] = 31;
    _atmo_chem.chempoint[40] = 32, _atmo_chem.chempoint[41] = 26;
    _atmo_chem.chempoint[42] = 5,  _atmo_chem.chempoint[43] = 6;
    _atmo_chem.chempoint[44] = 27, _atmo_chem.chempoint[45] = 12;
    _atmo_chem.chempoint[46] = 28, _atmo_chem.chempoint[47] = 30;
    _atmo_chem.chempoint[48] = 29, _atmo_chem.chempoint[49] = 7;
    _atmo_chem.chempoint[50] = 8,  _atmo_chem.chempoint[51] = 18;

    _atmo_chem.molar_mass[0]  = 30.000, _atmo_chem.molar_mass[1]  = 46.000;
    _atmo_chem.molar_mass[2]  = 48.000, _atmo_chem.molar_mass[3]  = 16.000;
    _atmo_chem.molar_mass[4]  = 16.000, _atmo_chem.molar_mass[5]  = 17.010;
    _atmo_chem.molar_mass[6]  = 33.010, _atmo_chem.molar_mass[7]  = 34.010;
    _atmo_chem.molar_mass[8]  = 62.010, _atmo_chem.molar_mass[9]  = 108.01;
    _atmo_chem.molar_mass[10] = 47.010, _atmo_chem.molar_mass[11] = 63.010;
    _atmo_chem.molar_mass[12] = 79.010, _atmo_chem.molar_mass[13] = 28.010;
    _atmo_chem.molar_mass[14] = 30.030, _atmo_chem.molar_mass[15] = 44.050;
    _atmo_chem.molar_mass[16] = 75.040, _atmo_chem.molar_mass[17] = 121.05;
    _atmo_chem.molar_mass[18] = 43.040, _atmo_chem.molar_mass[19] = 74.040;
    _atmo_chem.molar_mass[20] = 120.04, _atmo_chem.molar_mass[21] = 47.030;
    _atmo_chem.molar_mass[22] = 47.030, _atmo_chem.molar_mass[23] = 77.040;
    _atmo_chem.molar_mass[24] = 46.070, _atmo_chem.molar_mass[25] = 16.040;
    _atmo_chem.molar_mass[26] = 47.030, _atmo_chem.molar_mass[27] = 32.040;
    _atmo_chem.molar_mass[28] = 48.040, _atmo_chem.molar_mass[29] = 46.030;
    _atmo_chem.molar_mass[30] = 30.070, _atmo_chem.molar_mass[31] = 47.030;
    _atmo_chem.molar_mass[32] = 60.050, _atmo_chem.molar_mass[33] = 76.050;
    _atmo_chem.molar_mass[34] = 15.030, _atmo_chem.molar_mass[35] = 16.000;
    _atmo_chem.molar_mass[36] = 28.050, _atmo_chem.molar_mass[37] = 27.050;
    _atmo_chem.molar_mass[38] = 56.110, _atmo_chem.molar_mass[39] = 68.180;
    _atmo_chem.molar_mass[40] = 70.090, _atmo_chem.molar_mass[41] = 136.24;
    _atmo_chem.molar_mass[42] = 92.140, _atmo_chem.molar_mass[43] = 106.16;
    _atmo_chem.molar_mass[44] = 108.14, _atmo_chem.molar_mass[45] = 141.15;
    _atmo_chem.molar_mass[46] = 48.040, _atmo_chem.molar_mass[47] = 108.14;
    _atmo_chem.molar_mass[48] = 72.060, _atmo_chem.molar_mass[49] = 64.060;
    _atmo_chem.molar_mass[50] = 98.080, _atmo_chem.molar_mass[51] = 63.030;

    const char *label[] = {"NO",   "NO2",  "O3",    "O3P",
                           "O1D",  "OH",   "HO2",   "H2O2",
                           "NO3",  "N2O5", "HONO",  "HNO3",
                           "HNO4", "CO",   "HCHO" , "ALD2",
                           "C2O3", "PAN",  "ALDX",  "CXO3",
                           "PANX", "XO2",  "XO2N",  "NTR",
                           "ETOH", "CH4",  "MEO2",  "MEOH",
                           "MEPX", "FACD", "ETHA",  "ROOH",
                           "AACD", "PACD", "PAR",   "ROR",
                           "ETH",  "OLE",  "IOLE",  "ISOP",
                           "ISPD", "TERP", "TOL",   "XYL",
                           "CRES", "TO2",  "OPEN",   "CRO",
                           "MGLY", "SO2",  "H2SO4", "HCO3"};
    const char *name[]
      = {"species_no",   "species_no2",  "species_o3",    "species_o3p",
         "species_o1d",  "species_oh",   "species_ho2",   "species_h2o2",
         "species_no3",  "species_n2o5", "species_hono",  "species_hno3",
         "species_hno4", "species_co",   "species_hcho",  "species_ald2",
         "species_c2o3", "species_pan",  "species_aldx",  "species_cxo3",
         "species_panx", "species_xo2",  "species_xo2n",  "species_ntr",
         "species_etoh", "species_ch4",  "species_meo2",  "species_meoh",
         "species_mepx", "species_facd", "species_etha",  "species_rooh",
         "species_aacd", "species_pacd", "species_par",   "species_ror",
         "species_eth",  "species_ole",  "species_iole",  "species_isop",
         "species_ispd", "species_terp", "species_tol",   "species_xyl",
         "species_cres", "species_to2",  "species_open",  "species_cro",
         "species_mgly", "species_so2",  "species_h2so4", "species_hco3"};

    if (_atmo_chem.species_to_field_id == nullptr)
      CS_MALLOC(_atmo_chem.species_to_field_id, _atmo_chem.n_species, int);

    for (int ii = 0; ii < _atmo_chem.n_species; ii++) {
      int f_id = cs_variable_field_create(name[ii],
                                          label[ii],
                                          CS_MESH_LOCATION_CELLS,
                                          1);
      cs_field_t *f = cs_field_by_id(f_id);
      cs_add_model_field_indexes(f);
      _atmo_chem.species_to_field_id[ii] = f->id;
    }

  }
  //User defined chemistry using SPACK file and routines
  else if (_atmo_chem.model == 4)  {

    /* This function read the number of species, their molar mass
     * and creates variables */
    cs_atmo_declare_chem_from_spack();

    int  spack_n_species, n_photolysis;
    cs_ext_polyphemus_ssh_dimensions(&spack_n_species,
                                     &_atmo_chem.n_reactions,
                                     &n_photolysis);

    if (spack_n_species != _atmo_chem.n_species)
      cs_parameters_error
        (CS_ABORT_IMMEDIATE,
         _("In atmospheric chemistry inpu data"),
         _("The number of gaseous species read from the SPACK file\n"
           "is not equal to the one read in the SPACK source file."));
  }
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief This function set the file name of the SPACK file.
 *
 * \param[in] file_name  name of the file.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_chemistry_set_spack_file_name(const char  *file_name)
{
  if (file_name == nullptr) {
    _atmo_chem.model = 0;
    return;
  }

  _atmo_chem.model = 4;

  CS_MALLOC(_atmo_chem.spack_file_name,
            strlen(file_name) + 1,
            char);

  strcpy(_atmo_chem.spack_file_name, file_name);
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief This function sets the file name to initialize the aerosol library.
 *
 * \param[in] file_name  name of the file.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_chemistry_set_aerosol_file_name(const char  *file_name)
{
  if (file_name == nullptr) {
    _atmo_chem.aerosol_model = CS_ATMO_AEROSOL_OFF;
    return;
  }

  _atmo_chem.aerosol_model = CS_ATMO_AEROSOL_SSH;

  CS_MALLOC(_atmo_chem.aero_file_name,
            strlen(file_name) + 1,
            char);

  strcpy(_atmo_chem.aero_file_name, file_name);
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief This function declares additional transported variables for
 *        atmospheric module for the chemistry defined from SPACK.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_declare_chem_from_spack(void)
{
  assert(_atmo_chem.model == 4);

  if (_atmo_chem.spack_file_name == nullptr)
    bft_error(__FILE__, __LINE__, 0,
              _("Atmo chemistry from SPACK file: requires a SPACK file."));

  char line[512] = "";

  /* Open file */
  bft_printf("SPACK file for atmo chemistry:\n    %s \n",
             _atmo_chem.spack_file_name);

  FILE *file = fopen(_atmo_chem.spack_file_name, "rt");
  if (file == nullptr)
    bft_error(__FILE__, __LINE__, 0,
              _("Atmo chemistry from SPACK file: Could not open file."));

  int line_num = 0;

  /* Read "[species]" */
  while (true) {
    line_num++;
    if (fscanf(file, "%s\n", line) != 1)
      bft_error(__FILE__, __LINE__, 0,
                _("Atmo chemistry from SPACK file: Could not skip header."));

    if (strncmp("[species]", line, 512) == 0)
      break;
  }

  /* Read SPACK: first loop count the number of species */
  for (int i = 1; true; i++ ) {
    /* Read species */
    line_num++;
    if (fscanf(file, "%s\n", line) != 1)
      bft_error(__FILE__, __LINE__, 0,
                _("Atmo chemistry from SPACK file: Could not read line %d."),
                line_num);

    /* When reach [molecular_waight]: break */
    if (strncmp("[molecular_weight]", line, 512) == 0)
      break;
    else {
      _atmo_chem.n_species = i;
    }
  }

  /* Now allocate arrays */
  CS_MALLOC(_atmo_chem.species_to_field_id, _atmo_chem.n_species, int);
  CS_MALLOC(_atmo_chem.molar_mass, _atmo_chem.n_species, cs_real_t);
  CS_MALLOC(_atmo_chem.chempoint, _atmo_chem.n_species, int);

  /* Read SPACK: second loop Create variables and read molar mass */
  for (int i = 0; i < _atmo_chem.n_species; i++ ) {
    char name[512] = "";
    char label[512] = "";

    /* Read species name and molar mass */
    line_num++;
    if (fscanf(file, "%s %lf\n", line, &(_atmo_chem.molar_mass[i])) != 2)
      bft_error(__FILE__, __LINE__, 0,
                _("Atmospheric chemistry from SPACK file:\n"
                  "  could not read species name and molar mass, line %d."),
                line_num);

    /* The order is already ok */
    _atmo_chem.chempoint[i] = i+1; //FIXME ?

    /* Build name of the field:
     * species_name in lower case */
    strcpy(name, "species_");
    _strtolower(label, line);
    strcat(name, label);

    /* Field of dimension 1 */
    /* Give the original name as label */
    _atmo_chem.species_to_field_id[i]
      = cs_variable_field_create(name, line, CS_MESH_LOCATION_CELLS, 1);

    cs_field_t *f = cs_field_by_id(_atmo_chem.species_to_field_id[i]);
    cs_add_model_field_indexes(f);
  }
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief Deactivate chemistry initialization procedure
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_chemistry_initialization_deactivate(void)
{
  _init_atmo_chemistry = 0;
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief Check if the chemistry module needs initialization
 *
 * \return int value : 1 if needed, 0 if not
 */
/*----------------------------------------------------------------------------*/

int
cs_atmo_chemistry_need_initialization(void)
{
  return _init_atmo_chemistry;
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief Initialize chemistry cell kinetic rates arrays.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_chemistry_initialize_reacnum(void)
{
  const cs_lnum_t n_cells = cs_glob_mesh->n_cells;

  if (_atmo_chem.reacnum == nullptr)
    CS_MALLOC(_atmo_chem.reacnum,
              (cs_lnum_t)(_atmo_chem.n_reactions)*n_cells,
              cs_real_t);
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief Finalize atmospheric chemistry structures.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_chemistry_finalize(void)
{
  if (_atmo_chem.aerosol_model != CS_ATMO_AEROSOL_OFF)
    cs_atmo_aerosol_finalize();

  CS_FREE(_atmo_chem.reacnum);
  CS_FREE(_atmo_chem.species_to_field_id);
  CS_FREE(_atmo_chem.molar_mass);
  CS_FREE(_atmo_chem.chempoint);
  CS_FREE(_atmo_chem.conv_factor_jac);
  CS_FREE(_atmo_chem.spack_file_name);
  CS_FREE(_atmo_chem.aero_file_name);
  CS_FREE(_atmo_chem.chem_conc_file_name);
  CS_FREE(_atmo_chem.conc_profiles);
  CS_FREE(_atmo_chem.z_conc_profiles);
  CS_FREE(_atmo_chem.t_conc_profiles);
  CS_FREE(_atmo_chem.x_conc_profiles);
  CS_FREE(_atmo_chem.y_conc_profiles);
  CS_FREE(_atmo_chem.dlconc0);
  CS_FREE(_atmo_chem.species_profiles_to_field_id);
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief Print the atmospheric chemistry options to setup.log.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_chemistry_log_setup(void)
{
  if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] < 0)
    return;

  cs_log_printf(CS_LOG_SETUP,
                _("\n"
                  "Atmospheric chemistry options\n"
                  "-----------------------------\n\n"));

  if (cs_glob_atmo_chemistry->model == 0) {

    /* No atmospheric chemistry */
    cs_log_printf(CS_LOG_SETUP,
                  _("  No atmospheric chemistry\n\n"));

  }
  else if (   cs_glob_atmo_chemistry->model == 1
           || cs_glob_atmo_chemistry->model == 2
           || cs_glob_atmo_chemistry->model == 3) {

    /* Pre-defined schemes */
    cs_log_printf
      (CS_LOG_SETUP,
       _("  Atmospheric chemistry activated\n\n"
         "    Pre-defined scheme %12d\n\n"
         "      n_species: %18d (Number of species)\n"
         "      n_reactions: %16d (Number of reactions)\n"
         "      photolysis: %17s\n"
         "      frozen_gas_chem: %12s\n\n"),
       cs_glob_atmo_chemistry->model,
       cs_glob_atmo_chemistry->n_species,
       cs_glob_atmo_chemistry->n_reactions,
       cs_glob_atmo_chemistry->chemistry_with_photolysis ? "Yes": "No",
       cs_glob_atmo_chemistry->frozen_gas_chem ? "Yes": "No");

  }
  else if (cs_glob_atmo_chemistry->model == 4) {

    /* User defined SPACK chemistry */
    cs_log_printf
      (CS_LOG_SETUP,
       _("  Atmospheric chemistry activated\n\n"
         "    User-defined SPACK scheme\n\n"
         "      n_species: %18d (Number of species)\n"
         "      n_reactions: %16d (Number of reactions)\n"
         "      photolysis: %17s\n"
         "      frozen_gas_chem: %12s\n"
         "      Spack file: %17s\n"),
       cs_glob_atmo_chemistry->n_species,
       cs_glob_atmo_chemistry->n_reactions,
       cs_glob_atmo_chemistry->chemistry_with_photolysis ? "Yes": "No",
       cs_glob_atmo_chemistry->frozen_gas_chem ? "Yes": "No",
       cs_glob_atmo_chemistry->spack_file_name);

  }
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief Print the atmospheric aerosols options to setup.log.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_aerosol_log_setup(void)
{
  if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] < 0)
    return;

  cs_log_printf
    (CS_LOG_SETUP,
     _("\nAtmospheric aerosols options\n"
       "----------------------------\n\n"));

  cs_atmo_aerosol_type_t atm_aer_type = cs_glob_atmo_chemistry->aerosol_model;

  if (atm_aer_type == CS_ATMO_AEROSOL_OFF)
    cs_log_printf(CS_LOG_SETUP, _("  %s\n"),
                  cs_atmo_aerosol_type_name[atm_aer_type]);

  else if (atm_aer_type == CS_ATMO_AEROSOL_SSH) {

    /* Atmospheric chemistry with SSH */
    cs_log_printf
      (CS_LOG_SETUP,
       _("  Atmospheric aerosols activated\n\n"
         "    Global parameters\n\n"
         "      model: %22s (%s)\n"
         "      n_layer: %20d (Number of layers inside aerosols)\n"
         "      n_size:  %20d (Number of resolved aerosols)\n"
         "      init_gas_with_lib: %10s\n"
         "      init_aero_with_lib: %9s\n"
         "      namelist: %s\n\n"),
       cs_atmo_aerosol_type_enum_name[atm_aer_type],
       cs_atmo_aerosol_type_name[atm_aer_type],
       cs_glob_atmo_chemistry->n_layer,
       cs_glob_atmo_chemistry->n_size,
       cs_glob_atmo_chemistry->init_gas_with_lib ? "Yes": "No",
       cs_glob_atmo_chemistry->init_aero_with_lib ? "Yes": "No",
       cs_glob_atmo_chemistry->aero_file_name);
  }
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief initialize gaseous and particulate concentrations and aerosol number
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_chem_initialize_dlconc0(void)
{
  if (_atmo_chem.dlconc0 != nullptr)
    return;

  const int n_aer = _atmo_chem.n_size;
  const int nlayer_aer = _atmo_chem.n_layer;
  const int size = n_aer*(1+nlayer_aer);

  CS_MALLOC(_atmo_chem.dlconc0, size, cs_real_t);
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief initialize (allocate)  chemistry concentrations profiles
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_chemistry_initialize_conc_profiles(void)
{
  const int nespg   = _atmo_chem.n_species;
  const int _nbchmz = _atmo_chem.n_z_profiles;
  const int _nbchim = _atmo_chem.nt_step_profiles;
  const int size    = _atmo_chem.n_species*_nbchmz*_nbchim;

  if (_atmo_chem.conc_profiles == nullptr)
    CS_MALLOC(_atmo_chem.conc_profiles, size, cs_real_t);

  if (_atmo_chem.z_conc_profiles == nullptr)
    CS_MALLOC(_atmo_chem.z_conc_profiles, _nbchmz, cs_real_t);

  if (_atmo_chem.t_conc_profiles == nullptr)
    CS_MALLOC(_atmo_chem.t_conc_profiles, _nbchim, cs_real_t);

  if (_atmo_chem.x_conc_profiles == nullptr)
    CS_MALLOC(_atmo_chem.x_conc_profiles, _nbchim, cs_real_t);

  if (_atmo_chem.y_conc_profiles == nullptr)
    CS_MALLOC(_atmo_chem.y_conc_profiles, _nbchim, cs_real_t);

  if (_atmo_chem.conv_factor_jac == nullptr)
    CS_MALLOC(_atmo_chem.conv_factor_jac, nespg*nespg, cs_real_t);
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief  Computes the explicit chemical source term for atmospheric
 *         chemistry in case of a semi-coupled resolution
 *
 * \param[in]     f_id          field id
 * \param[out]    st_exp        explicit part of the source term
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_chem_exp_source_terms(int          f_id,
                              cs_real_t    st_exp[])
{
  const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
  const cs_real_t *cell_vol = cs_glob_mesh_quantities->cell_vol;

  const cs_atmo_chemistry_t *atmo_chem = cs_glob_atmo_chemistry;
  const int n_reactions = atmo_chem->n_reactions;
  const int n_species_g = atmo_chem->n_species;

  const cs_real_t navo = cs_physical_constants_avogadro;

  /* If gaseous chemistry is computed with the external library SSH-aerosol
   * The source term was estimated at the beginning of
   * cs_solve_transported_variables.
   * We update the source term and return directly */

  if (atmo_chem->aerosol_model != CS_ATMO_AEROSOL_OFF)
    cs_parameters_error
      (CS_ABORT_IMMEDIATE,
       _("Partially coupled chemistry combined with external"),
       _(" aerosol library not implemented yet\n"));

  cs_real_t conv_factor[n_species_g], source[n_species_g];
  cs_real_t rk[n_reactions], dlconc[n_species_g], dchema[n_species_g];

  const cs_real_t *cpro_rho = CS_F_(rho)->val;
  const cs_real_t *reacnum = atmo_chem->reacnum;

  cs_real_t **cvara_espg = nullptr;
  CS_MALLOC(cvara_espg, n_species_g, cs_real_t *);

  /* Arrays of pointers containing the fields values for each species */
  for (int spe_id = 0; spe_id < n_species_g; spe_id++) {
    cs_field_t *f = cs_field_by_id(atmo_chem->species_to_field_id[spe_id]);
    cvara_espg[spe_id] = f->val_pre;
  }

  const int idx_chema
    = atmo_chem->chempoint[f_id - atmo_chem->species_to_field_id[0]] - 1;

  for (int c_id = 0; c_id < n_cells; c_id++) {

    const cs_real_t rho = cpro_rho[c_id];

    // Filling working array rk
    for (int ii = 0; ii < n_reactions; ii++)
      rk[ii] = reacnum[ii*n_cells+ c_id];

    // Filling working arrays
    for (int spe_id = 0; spe_id < n_species_g; spe_id++) {
      const int idx = atmo_chem->chempoint[spe_id] - 1;
      source[spe_id] = 0.0;
      dlconc[idx] = cvara_espg[spe_id][c_id];
      conv_factor[idx] = rho*navo*(1.0e-9)/atmo_chem->molar_mass[spe_id];
    }

    // Computation of C(Xn)
    if (atmo_chem->model == 1)
      cs_ext_polyphemus_fexchem_1(n_species_g,
                                  n_reactions,
                                  dlconc, rk, source,
                                  conv_factor, dchema);
    else if (atmo_chem->model == 2)
      cs_ext_polyphemus_fexchem_2(n_species_g,
                                  n_reactions,
                                  dlconc, rk, source,
                                  conv_factor, dchema);
    else if (atmo_chem->model == 3)
      cs_ext_polyphemus_fexchem_3(n_species_g,
                                  n_reactions,
                                  dlconc, rk, source,
                                  conv_factor, dchema);
    else if (atmo_chem->model == 4)
      cs_ext_polyphemus_fexchem_4(n_species_g,
                                  n_reactions,
                                  dlconc, rk, source,
                                  conv_factor, dchema);

    /*Adding source term to st_exp
     * The first n_species_g user scalars are supposed to be chemical species
     * TODO: try to implicit the ST */
    st_exp[c_id] += rho*cell_vol[c_id]*dchema[idx_chema];
  }

  CS_FREE(cvara_espg);
}

/*----------------------------------------------------------------------------*/
/*!
 * \brief  Calls the rosenbrock resolution for atmospheric chemistry
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_compute_gaseous_chemistry(void)
{
  const cs_lnum_t n_cells = cs_glob_mesh->n_cells;

  const cs_atmo_chemistry_t *atmo_chem = cs_glob_atmo_chemistry;
  const int n_reactions = atmo_chem->n_reactions;
  const int n_species_g = atmo_chem->n_species;
  const cs_real_t dtchemmax = atmo_chem->dt_chemistry_max;

  const cs_real_t navo = cs_physical_constants_avogadro;

  const cs_real_t *dt = CS_F_(dt)->val;
  const cs_real_t *cpro_rho = CS_F_(rho)->val;
  const cs_real_t *reacnum = atmo_chem->reacnum;

  cs_real_t conv_factor[n_species_g], source[n_species_g];
  cs_real_t rk[n_reactions], dlconc[n_species_g], dchema[n_species_g];

  cs_real_t **cvar_espg = nullptr;
  cs_real_t **cvara_espg = nullptr;
  CS_MALLOC(cvar_espg, n_species_g, cs_real_t *);
  CS_MALLOC(cvara_espg, n_species_g, cs_real_t *);

  /* Arrays of pointers containing the fields values for each species */
  for (int spe_id = 0; spe_id < n_species_g; spe_id++) {
    cs_field_t *f = cs_field_by_id(atmo_chem->species_to_field_id[spe_id]);
    cvar_espg[spe_id] = f->val;
    cvara_espg[spe_id] = f->val_pre;
  }

  for (int c_id = 0; c_id < n_cells; c_id++) {

    const cs_real_t dtc = dt[c_id];
    const cs_real_t rho = cpro_rho[c_id];

    // Filling working array rk
    for (int ii = 0; ii < n_reactions; ii++)
      rk[ii] = reacnum[ii*n_cells+ c_id];

    // Filling working arrays
    for (int spe_id = 0; spe_id < n_species_g; spe_id++) {
      const int idx = atmo_chem->chempoint[spe_id] - 1;
      source[spe_id] = 0.0;
      conv_factor[idx] = rho*navo*(1.0e-9)/atmo_chem->molar_mass[spe_id];
    }

    if (   atmo_chem->chemistry_sep_mode == 1
        || cs_glob_time_step->nt_cur < cs_glob_time_step->nt_ini) {
      /* split Rosenbrock solver
         ----------------------- */

      // Filling working array dlconc with values at current time step
      for (int spe_id = 0; spe_id < n_species_g; spe_id++) {
        const int idx = atmo_chem->chempoint[spe_id] - 1;
        dlconc[idx] = cvar_espg[spe_id][c_id];
      }

    }
    else {
      /* semi-coupled Rosenbrock solver
         ------------------------------  */

      for (int spe_id = 0; spe_id < n_species_g; spe_id++) {
        const int idx = atmo_chem->chempoint[spe_id] - 1;
        dlconc[idx] = cvara_espg[spe_id][c_id];
      }

      // Computation of C(Xn)
      if (atmo_chem->model == 1)
        cs_ext_polyphemus_fexchem_1(n_species_g,
                                    n_reactions,
                                    dlconc, rk, source,
                                    conv_factor, dchema);
      else if (atmo_chem->model == 2)
        cs_ext_polyphemus_fexchem_2(n_species_g,
                                    n_reactions,
                                    dlconc, rk, source,
                                    conv_factor, dchema);
      else if (atmo_chem->model == 3)
        cs_ext_polyphemus_fexchem_3(n_species_g,
                                    n_reactions,
                                    dlconc, rk, source,
                                    conv_factor, dchema);
      else if (atmo_chem->model == 4)
        cs_ext_polyphemus_fexchem_4(n_species_g,
                                    n_reactions,
                                    dlconc, rk, source,
                                    conv_factor, dchema);

      /* Explicit contribution from dynamics as a source term:
         (X*-Xn)/dt(dynamics) - C(Xn). See cs_polyphemus_spack.f90.
         The first n_species_g user scalars are expected to be
         chemical species */
      for (int spe_id = 0; spe_id < n_species_g; spe_id++) {
        const int idx = atmo_chem->chempoint[spe_id] - 1;
        source[idx] = (cvar_espg[spe_id][c_id] - cvara_espg[spe_id][c_id])/dtc
                    -  dchema[idx];
      }
    }

    // Rosenbrock resolution

    // The maximum time step used for chemistry resolution is dtchemmax
    if (dtc < dtchemmax) {
      _rosenbrock_solver(dlconc, source, source, conv_factor, dtc, rk, rk);
    }
    else {
      const int ncycle = (int)(dtc/dtchemmax);
      const int rest = (int)dtc%(int)dtchemmax;
      const cs_real_t dtrest = rest;
      for (int ii = 0; ii < ncycle; ii++)
        _rosenbrock_solver(dlconc,
                           source,
                           source,
                           conv_factor,
                           dtchemmax, rk, rk);

      _rosenbrock_solver(dlconc,
                         source,
                         source,
                         conv_factor,
                         dtrest, rk, rk);
    }

    // Update of values at current time step
    for (int spe_id = 0; spe_id < n_species_g; spe_id++) {
      const int idx = atmo_chem->chempoint[spe_id] - 1;
      cvar_espg[spe_id][c_id] = dlconc[idx];
    }

  } // end loop on cells

  CS_FREE(cvar_espg);
  CS_FREE(cvara_espg);

  // Clipping
  for (int spe_id = 0; spe_id < n_species_g; spe_id++) {
    cs_field_t *f = cs_field_by_id(atmo_chem->species_to_field_id[spe_id]);
    cs_scalar_clipping(f);
  }
}

/*----------------------------------------------------------------------------*/

END_C_DECLS
